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Abstract 



In this paper we introduce efficient Monte Carlo estimators for the val- 
uation of high- dimensional derivatives and their sensitivities ("Greeks"). 

QQ , These estimators are based on an analytical, usually approximative rep- 

resentation of the underlying density. We study approximative densities 
obtained by the WKB method. The results are applied in the context of 

P^ . a Libor market model. 
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1 Introduction 



Valuation methods for high-dimensional derivative products are typically based 
on Monte Carlo simulation of the underlying process. The dynamics of the 
^N^ , underlyings are usually given via a (jump-) diffusion SDE. In case of a diffu- 

sion SDE, the underlying process may be simulated using an Euler scheme or 
r — ^ , a (weak) second order scheme e.g. see Kloeden & Platen [23| or Milstein & 

(^ ' Tretyakov [29|. For simulation of jump-diffusions see e.g. Cont & Tankov [7|, 

OO , and Glasserman & Merener [l6| for simulation of (Libor) interest rate models 

^^ ' with jumps. 

The evaluation of option sensitivities, 'Greeks' in financial terms, comes 
ij - down to the computation of expressions of the form -^E{f{X^)) (and possibly 

rS , higher order derivatives), where / is a pay-off function, X is the state of an 

d ' underlying process depending on some parameter A. For example, the first and 

second order derivatives with respect the initial state are called Deltas and Gam- 
mas, respectively. In the literature the evaluation of Greeks has been treated 
by several methods (a nice overview about classical and recent literature is pro- 
vided in Elie, Fermanian, Touzi [Uj). Classical finite difference approaches have 



been studied by L' Ecuyer & Perron [10|, Broadie & Glasserman [4], Milstein & 



ypro 
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Schoenmakers [28|, Milstein & Tretyakov [30|, Detemple, Garcia, Rindisbacher 



I , and Giles & Glasserman [l5| . These approaches are quite general and easy 
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to implement as they do not require particular knowledge of the distribution of 
the underlying. However, they require full blown simulation of the correspond- 
ing system of stochastic differential equations and, in order to be efficient, some 
degree of regularity with respect to the pay-off function. In case the transition 
kernel of X is known or known in a good approximation, the latter drawback can 
be avoided by differentiating this kernel with respect to the sensitivity param- 



eter A, see Fries & Kampen [1J| and Fries & Joshi 13]. The typical difficulty 
that the distribution of the underlying is only known for very special cases 
was overcome by Fournie, Lasry, Lebuchoux, Lions, Touzi [12| . who used the 
Malliavin integration-by-parts formula in order to express Greeks in the form 
E{f{X)Tr), where the random variable tt is called a Greek weight. In a more 
recent alternative approach Elie, Fermanian, Touzi [llj construct Greek esti- 
mators which are based on variance minimizing choices of Greek weights. As 
a matter of fact, the Malliavin method does not lead to this optimal weight in 
general. In order to avoid straightforward SDE simulation in the context of the 
Libor market interest rate model, and so reducing simulation costs, Kurban- 
niuradov, Sabelfeld & Schoenmakers [26| considered lognormal approximations 
for the transition density, whereas Hunter, Jackel, & Joshi [12], and Pelsser, 
Pietersz & van Regenmortel [3l| propose specific drift approximations. 

In an ideal situation, the density of the underlying process X^ at a fixed 
point in time is known explictly and an efficient method to sample from it 
is available. Usually, however, neither of this is true. Even if the transition 
density is known, we will show that calculation of sensitivities, based on kernel 
differentiation for instance (as in ,14j). may cause problems (high variance) in 
case the kernel under consideration is 'highly peaked ', for example due to small 
maturities, low volatilities, or high dimensionality of the underlying system. In 
this paper we therefore choose for a rather general approach with the following 
objectives. 

• Developing efficient variance bounded probabilistic representations for price 
sensitivities, based on an analytical approximation of the underlying den- 
sity and a possibly rougher approximative standard density (e.g. a log- 
normal density) which is basically used as an importance sampler. 



• 



Construction of a " good" analytical approximation for the density of the 
underlying process by using (convergent) WKI^^ methods; 



We underline that, in principle, the way of constructing an analytical approx- 
imation of the transition density is not essential for the developed Greek esti- 
mators. In this article we exploit the use of WKB approximations as a generic 
convergent method. In special cases, however, construction of high accuracy 
transition kernels may be possible by other means (see [26|] for example). 

The structure of the paper is as follows. In Section 2 we set up the model 
class for which we exemplify our methods and specify the financial products 
(including Bermudan callables) for which prices and sensitivities are to be de- 
termined. In Section 3 we introduce probabilistic representations for integral 
functionals of kernel type and their derivatives. As a particular result we prove 
that the corresponding estimator for the derivatives has non-exploding variance 



^The historical origin of the name is the work of Wentzel, Kramers, BriouUin in the 
context of semiclassical solutions of the Schrodinger equation. The meaning of WKB has 
broadened since; nowadays, it refers to analytic expansions of exponential form. 



for sharply peaked kernels in contrast to some existing weighted Monte Carlo 
schemes. This estimator thus allows for efficient Monte Carlo estimation of 
option sensitivities, in particular with respect to under lyings (Deltas), even in 
situations where the densities are sharply peaked (for instance when volatilities 
are small). The general probabilistic representations introduced in Section 3 are 
applied to the computation of Deltas for Bermudan callable products in Sec- 
tion 4. Section 5 deals with the WKB-theory of densities of diffusion equations 
(densities of processes which have continuous paths). In Section 5.1 we summa- 
rize some results concerning pointwise valid WKB-representations of densities 



obtained in Kampen [2]J . Since in practice only finitely many terms of a WKB 
expansion can be computed, it will be necessary to use a truncated form of 
the WKB-representation for actual computations. In Section 5.2. we analyze 
the effect of this truncation error on approximations of solutions of Cauchy 
problems and their derivatives. The case of non-autonomous diffusion models is 
discussed in Section 5.3. The results of Sections 2-5 are applied in Section 6 to 
the Libor market model. In Section 6.1. we compute explicitly the first three 
coefficients of the WKB representation of the Libor model density. In Section 
6.2 we compute prices and Deltas in a case study of European swaptions. 

2 Basic setup 

Let X ~ {X^, ..., X") be a Markovian process of financial derivative in M" (]R_|_ 
:= {x : X > 0}) under a given pricing measure P, connected with a given dis- 
counting numeraire B, B > 0, on some filtered probability space. For example, 
X may represent a system of asset prices or (Libor) interest rates. A popular 
framework for the system (X, B) is, for instance, the class of jump-diffusions 
(e.g. Cont & Tankov 0)- For simplicity however, we mainly consider in the 
present article ordinary diffusions, but, note that the main results generally 
extend to jump processes as well (see Kampen, Kolodko, Schoenmakers |22l|). 
With respect to an n-dimensional standard Wiener process W = {W^, ..., 
M^")^ on the probability space {fl,!F, i^t)te[to,T]TP): where as usual (J^t) is the 
P- augment at ion of the filtration generated by W, we assume that X is governed 
by the stochastic differential equation (SDE), 

j-^f,it,X)dt + Y,'^'Ht,X)dW^, l<i,3<n. (1) 

3=1 

It is assumed that /x(i, x) and the matrix a{t, x) = (c'-' (t, x)) , t G [fg, T], x € K" 
are such that for all xq € M" , there exists a unique solution t ^ Xt € K" of 
dl]) for to < i < T satisfying Xf,, = xo =: Xl°'^° . It is further assumed that the 
Markov process X has a transition density 

p(t,a;,s,y), h<s<t<T, x,yeR!^, (2) 

which is differentiable with respect to x, y, s, and t, up to any order. In order to 
guarantee the existence and uniqueness of ([T]), and the existence of the transition 
density Q as stated, it is sufficient to require that the functions /x(-, •) and 
(t(-, •) are bounded and have bounded derivatives up to any order, and that the 
volatility matrix cr(i, x) is regular with 

< Ai < l(crcr^) {t,x)\ < A2 (3) 



for all (t,x), t E [to,T], x E M", and some < Ai < A2 (see for example Bally 
& Talay [3]). 

Let us take (w.l.o.g.) Bq = 1 and consider a contingent claim with pay- 
off function of the form /{X^-jBr at some ( JF. )-stopping time t. By general 
arguments (e.g. Duffie [9(), the price of this claim at time to is given by 

u{to,xo)^E fiXl«^^«). 

For deterministic t, say t = T, we have a European claim, and for tg < t < T 
its discounted value process can be represented by 

ut := uit,Xt) :- E^^fiXr) = J p{t,Xt,T,y)f{y)dy, where 



u{t,x) = / p{t,x,T,y)f{y)dy (4) 

is the unique solution of the Cauchy problem 

i,j — 1 i— 1 

wCTja;) = f{x). 

The density kernelp(-, •, T, y) is the unique (weak) solution of ([5|) withp(T, a;, T, y) 
= 5(x — y), where 5 is the Dirac-delta function in Schwarz distribution sense. 

Of particular importance are Bcrmudan callable contracts. A Bermudan 
contract starting at io, is specified by a set of exercise dates {^1,^2, •■•,^1}, 
where to < ti < ... < ti <T, and corresponding (discounted) pay-off functions 
fi{x), 1 < i < I. According to the contract, the holder has the right to call 
(once) a cash-flow fi{X^°'^°)Bf°'^°' (with Bq°'^°' = 1) at an exercise date ti 
of his choice. It is well known (e.g. Duffie |9|]) that the discounted price of this 
contract at time t, to <t <T, assuming that no exercise took place before t, is 
given by 

uit,x) := sup £;/r(X*'") = Ef^t.4X%), t,_i < t < t„ (6) 

where x = X^"'^" , Tl^j the set of stopping times r taking values in {ti, t^+i, ..., tj}, 
and r»'^ is an optimal stopping time. In particular the process u{t,Xt) is a su- 
permartingale and is called the Snell envelope of the (discounted) cash-flow 
process Si{Xt.). 

3 Probabilistic representations and their esti- 
mators 

In this section we consider for a given smooth function u : M" -^ K+ and a 
smooth kernel function p : R" x R" -^ R+ , probabilistic representations for the 
integral 

I{x) :— I p{x, y)u{y)dy, and its gradient 
^(-) = J Q^p{-,yHy)dy, with _:=(^_,...,_j. (7) 



Here and in the following sufficient (uniform) integrability conditions are as- 
sumed to be fulfilled, for instance, in order to guarantee that ([7]) is valid. 

Remark 1. In ([7]), kernel p (which may or may not be a density in the second 
argument) and function u have to be distinguished from the respective defini- 
tions in Section m although they may be related. For fixed t,T, < t < T, one 
could take (see (|4][5])), p{x,y) := p{t,x,T,y) and u{x) := u{t,x) for example. 

Let C be an R"-valued random variable on some probability space with 
density </),</)> 0. Then, obviously, 

I(^)^Ep{x,Oj^ (8) 

is a probabilistic representation for ([7]) which may be estimated by the unbiased 
Monte Carlo estimator 

where for m = 1, ...,M, mC ^-^e i.i.d. samples from a distribution with density 
(j). By taking gradients in ([5]) we readily obtain the probabilistic representation 

with corresponding estimator, 

dl 1 s-^ d ..uirnO .-.-.x 

m—l ^ ' 

While as a rule Q is an effective estimator for I{x) for a proper choice of 
(p, unfortunately the gradient estimator ()lip has a serious drawback: If the 
kernel p{x,-) is sharply peaked (nearly proportional to a 'delta- function'), its 
variance may be extremely high. This fact is demonstrated by the following 
stylistic example of a multi-asset model, which is nevertheless realistic in orders 
of magnitude. 

Example 2. Consider for fixed xq G M", parameters s > 0, and ct > 0, the 
n-dimensional lognormal density 



2 " exp 

p(s,a;xo,y) := ^^_ ^ ^^/^ II 



-O^ln^iC 

2a '^ s Xn 



(27ra2s)"/'- 



(12) 



In p^ p{s, cr; xo, •) is the density of the random variable {x\e'^^^ , ..., x^e'^^^"), 
where C i i — l,...,rf, are i.i.d. standard normal random variables. Thus, for 
small s and a, p{s,a;xo, ■) is peaked ('delta-shaped') around xq. Let us now 
take (/)(•) := p{s,a;xo,-) in ([H]) and (fTU|) . respectively, and m = ||xo|| (a constant 
of order xq in magnitude). Clearly, estimator ([9]) equals ||xo|| almost surely and 



so has zero variance. However, estimator (jlip is not deterministic and we have 

M 



T^i'^o) := X7 E 



Ikol 



d 



M ^^ p{s, a; xo,m, C) dx^ 
INoll Y^ 5 , / /-N 

m—1 



■p{s,a;xo,mO 



\xo\_ 
M 



M In. 



,C^ 



E 



M 



M 

E 



.e^ 



„^, ^252;)^ ivi —^ a^/sx: 



Hence, E 



di 

dxi 



(xo) 



= as should be, but. 



Var 



dl 



(xo) 



\^o/4\\' 1 

M ct2 



(13) 



which explodes when a'^s goes to zero! 

Remark 3. In Fries & Kampen [1J| estimators (O and ([TT|) are used for com- 
puting prices and sensitivities of European Libor options, respectively. In their 
numerical examples they used 50% (rather high) volatility in order to amplify 
Monte Carlo errors. While, indeed, a larger volatility generally gives rise to a 
large Monte Carlo error of ^, Example [H shows that the opposite is true for es- 
timator (jlip . For example, 50% volatility in combination with 0.5 yr. maturity 
corresponds to a (just moderate) variance factor 1/ (cr^s) = 8.0 in (fT^. while 
a more usual Libor volatility, e.g. 14%, and 0.5 y maturity would give a factor 
102.0(!). 

In the present paper we propose sensitivity estimators which are efficient on 
a broad time and volatility scale. As a result, the next theorem provides a tool 
for constructing sensitivity (gradient) estimators with non-exploding variance. 

Theorem 4. Let X be a reference density on R" with X{z) ^ for all z (for 
example, the standard normal density). Let ^ be an W^ -valued random variable 



on some probability space, with density A and g : M" x 



be a smooth 



enough map which has at least continuous derivatives with \dg{x, z)/dz\ 7^ 0, 
such that for each x £ R" the random variable C^ := g{x,S,) has a density 4>{x, ■) 
on R" . Then, for ^ we have the probabilistic representation 



dL 



{x)=E 



d pix,CMC) 



E 



d p{x,g{x,£,))u{g{x,0) 



dx dx (j){xX^) dx 

with corresponding Monte Carlo estimator 



(t>ix,g{x,^)) 



— ix) = J- V ^ Pi^^9ix,„iO)u{g{x,mS.)) 



(14) 



(15) 



-^^^ II'IIq •— \/i? I 'I" where 
matrix norm. Then it holds 



denotes either a vector norm or a compatible 



E 



d p{x,g{x,£,))u{g{x,£,)) 



dx 



(t>{x,g{x,S,)) 



< 2mImImI 



- amImImI 



amImImImI, 

(16) 



hence the second moments of the Monte Carlo samplers for the components of 
81/ dx are bounded by the right-hand- side of 17^) . if for fixed x e R" , there are 
constants ai, ..., ag > 1 and Mi, ..., Mg > with 

111 111 1111 

h 1 = 1, \ \ = 1, \ 1 \ = 1, 

a4 ai a^ a^ a2 a^ a^ ai ae as 



such that, 



\\u{gix,m\2a,<Mi, 
99, 



I dp Idcj) 
p dx <j) dx 






<Mn 



p{x,g{x,£,)) 



<Mo 



2q2 



<Mr 



4>{x,g{x,C)) 



\dp 1 



< Ma 



2q4 



p dy 4> dy 



{x,g{x,^)) 



(17) 



<Mfi 



2q:6 



Proof. For any bounded measurable -0 : M" ^ R, we have 

dg{x,z) 



Tjj{g{x,z))X{z)dz = / %lj{g{x,z))(j){x,g{x,z)) 



dz 



dz. 



Therefore, the densities 4> and g are connected via the relationship 

dg{x,z) 



(j){x,g{x,z)) 



dz 



X{z) 



(18) 



By dUl), the right-hand-side of UH) equals 

d T:.p{x,9{x,S,))u{g{x,£,)) _ 9 f p{x,g{x,z))u{g{x,z)) 



- W 

dx (j){x,g{x,£,)) 



dx 
d_ 
dx 



(f){x,g{x,z)) 
p{x,g{x,z))u{g{x,z)) 



X{z)dz 



dg{x,z) 



dz 



dz 



To prove the moment estimation (J16p . we observe that 



E 



d pix,CHC) 



dx (j){xX^) 
E 



= E 



d p{x,g{x,£,))u{g{x,£,)) 



dx 4>{x,g{x,^)) 



I ( cY. 9 p{x,g{x,^)) p{x,g{x,^)) du ^..dg 

dx(p[x,g[x,i)) (p{x,g{x,Q) dy dx 



< 2E 



p^{x,g{x,0) 

<f>'^ix,g{x,£,)) 

+ 2Eu\g{x,0) 



i9{x,0) 



du 
dy 

d p{x,g{x,£,)) 



dx 

2 



{^,0 



dx(l){x,g{x,C)) 
Then by Holders inequality, (/) < 



= :2(/) + 2(//). 



IE 



du 
dy 



(5(2^,0) 



2(22 



'E 



dg_ 
dx 



(^,0 



2Q3 



J^ p^<-^{x,g{x,i)) 



IE 



(/.2"4(a;,5(x,C)) 



< Mi Mi Mi. 



For the second term we have (//) = 



Eu\g{x,0) 



p^{x,g{x,£,)) 



4>'^{x,g(x,S,)) 

%{x,g{x,0) if(x,g(x,0) 
p{x,g{x,S,)) (j)(x,g{x,£,)) 



dp 



^ix,g{x,0) ^{x,gix,0) 

Hx,9{x,£,)) 
1 



pix,gix,0) 



dg , ^^ 



<2Eu\g{x,0) 



2Eu^g{x,0) 



P^{x,g{x,0) 
<l>'^{x,g{x,£,)) 

p'^{x,g{x,i)) 



||(x,5(x,a) i!(^,5(x,0) 



p{x,g{x,S,)) 4>ix,g{x,^)) 
^{x,gix,0) ^(.x,gix,0) 



p{x,g{x,C)) <t^{x,g{x,i)) 



^9 ( c\ 



(l)'^{x,g{x,0) 
< 2MlMlMl + 2MlMlMlMl, 

again by Holders inequaUty. ■ 

Remark 5. If in PT)) the random variables M(.g(a;,^)), ^(^(x,^)), and so on, 
have moments of high enough order, Theorem [H guarantees that the variance of 
estimator (fTSJ) is controlled via the moment estimates (|17p . The most delicate 
bound in ([TT]) is A/5 in fact. Indeed, if one takes g{x,£,) = 5(2:0, estimator 
(fT5)) coUapses to ([TT|) . and in Example [U page [3 where 4>{x,y) = p{xo,y) in 
fact, we see that M^ cannot be taken small when a^s is small, i.e. when p 
is highly peaked around x. In contrast, if for fixed x, 4>{x^-) is approximately 
proportional iop{x, •) and 91n0(x, ■)/dx ~ d\Tip[x, ■)/dx (both with respect to 
the weight function (f){x, •)), a small M5 may exist. Note that for 0(-, •) exactly 
proportional to p{-, •), we may take M5 = 0. 

Remark 6. It can be shown that Theorem [H can be extended to probabilistic 
representations and corresponding estimators for higher order derivatives, 



dl 



{x)=E 



d p{x,Q^)u{C,^ 



E 



d p{x,g{x,C))u{g{x,^)) 



dxP^ ' dxP (/)(a;,C^) dx^ 

with corresponding Monte Carlo estimator 



Hx,g{x,0) 



M 



^1 .± 

dxP ^^^ M ^ dxP 

m— 1 



d p{x,g{x,^£,))u{g{x,„iO) 



4>{x,g{x,rnO) 



(19) 



where /3 := (/3i, . . . , /?„), (3i £ {0, 1,2,.. .} is a multi-index with (formally) dx^ — 
dx{^dx2 ■ ■ ■ 9x^" . Loosely speaking, the variance of the higher order derivative 
estimator (|19p can be bounded from above by an expression like (|16p involving 
(i) sufficiently high moments of the derivatives, y — > ^^, and 



dg(x,z) 



for 



fixed x, 7 < /3 (component wise), with respect to weight functions y — > <p{x, y) 
and z -^ A(z), respectively, and, (ii) for fixed x, L'^(IR" , 0(x, y)(i2/)-norms of 



d 

dx^ 



(f'{x,y) 
p{x,y) 



and 



d ((j){x,y) 



dy-y \p{x,y) 



7</3, 



for q large enough. 



Remark 7. In the next section we consider financial applications where I{x) 
is the price of a derivative contract considered in dependence of the argument x 
which may stand for the underlying process or some parameter (vector) which 
affects the dynamics of the underlying process (e.g. volatilities). Moreover, 
we there give a recipe how to construct a lognormal density approximation (/>, 
corresponding to a particular normal reference density and an exponential type 
transformation (in Theorem 31 A and g respectively). 

Remark 8. Theorem 4 can be can be generalized to the case where both p 
and (f) live on a common submanifold rather than on the whole state space. 
Via such a generalization it would be possible to extend our results to factor 



reduced situations in the spirit of Fries & Kampen [ij] and Fries & Joshi [13| . 
However, this is considered beyond the scope of the present article and therefore 
we restrict ourselves to the full-factor case. 



4 Sensitivities for Bermudan options 

Theorem H] may be applied in general for computing sensitivities (" Greeks" ) 
of derivative products. For estimator ()11|1 the danger of exploding variance 
is typically the largest when derivatives of prices with respect to underlyings 
(Deltas, Gammas) are considered. We therefore consider in this section only 
(first order) derivatives with respect to the underlying process, hence Deltas. 

Let T : n ^ M+ be a given stopping time with respect to the filtration (JF. ) . 
As usual we may think of fl as being the space of functions oj : [0, cxd) -^ M", 
which are continuous from the right and have limits from the left, and define 
T^'^{uj) := s + t{X^'^,Juj)). We now consider the Bermudan contract introduced 

in Section H For fixed t,t+ , to < t < t+ < h, x e W^, we have r*''^ = t* ' '"^ 
since t^'^ > ti, and we thus may write 



u{t,x) := EfiXlZ) - EE^^+fix]l%) 



p{t, X, t+,y)dyEfiX'^^fJ = j p{t, x, t+ ,y)u{t+ ,y)dy, 

by the Chapman-Kolmogorov equation. 

For each t, i+ as above, let 0(t, x, i+, y), g{t,x,t~^,y), and reference density 
X{t, i+, z) be as in Theorem 21 We then have the probabilistic representation 

^ ' <P{t,x,t+,g{t,x,t+,Oy^ ^*+,«(t,...+ .o^ V ; 

with Monte Carlo estimator 

JT/'y- ^^ — "*" V^ P{iiX,t+ ,g{t,X,t+ ,raS,)) ^ ^ ^1+ .g{t,x,t+ . „C) ^ /oi \ 

- M^^<^(t,:.,iH-,g(i,x,t+,™0)^^V-<*'-+' -'^' ^''^ 
and for the gradients (Deltas) we have the probabilistic representation 

A —^"r^^^-p ^ / p(t, X, t+,g{t, X, t+, 0) ^t+Mt.x.t+.i) A , . 



with Monte Carlo estimator 

A _ 1 Y^ g f p{t,x,t+,g{t,x,t+,„,^)) t+,git,^,t+,^i).\ ,„„^ 

where m^, to = 1, ...,M, are i.i.d. samples from the reference density A. Indeed, 
by pre-conditioning on J^j+ and then taking expectations we see that (PTj) and 
(f23|) are unbiased Monte Carlo estimators for the price (|20p and 'deltas' ([22|) . 
respectively. Moreover, if is close to p in the sense of Theorem |4j it is not 
difficult to see that also gradient estimator (P5|) has non-exploding variance 
when i+ J, t. 

Estimators ()2ip and (j23p are useful if one has an analytic approximation 
p{t, X, i+, y) of the density p(i, x, t'^,y) and known densities (l){x, •) for x S R" . 
The approximation p may be obtained by some specific method, for example by 
a WKB expansion as presented in Section [51 or some lognormal approximation 
as proposed in Kurbanmuradov, Sabelfeld & Schoenmakers [26[ for the Libor 
market model. Of course the density has to be chosen with some care. If it is 
possible to sample directly from p (e.g. in case of a log-normal approximation) 
we may take (j) = p. li not, (e.g. in the case of a WKB expansion) one may take 
for (j) a (not necessarily very accurate) lognormal approximation of the density 

P- 

A canonical lognormal approximation for p{t, x, t^ , z) is obtained by freezing 
X in the coefficients of ^ at the initial time. We thus obtain 

X^'°*^^'* :=a;'exp --^ / {G'^f{s,x)ds^ / r{s,x)ds 

+ 11/ ^'^(s,x)dVFi =:a;,exp(e^). (24) 

Here, (Ci)"=i is a Gaussian random vector with 

Eii = — V/ {a'^f{s,x)ds+ r(s,a;)ds=:/i^'*'*'''^, 1 < i < n, 

2 ^^, Jt Jt 

and 

Cov{^,,Q=Y. <j'\s,x)<j^'{s,x)ds^:a'r^'''*'^ l<i,j<n. 

Clearly, the density (j) is then given by 

12 Tl 

fj^ 1 1+ x 1 1+ X (In ^, In -T-, ..., In ■^) 

y* > 0, 1 < i < 71, with ^ t,t+ ,x t,t+.x being the density of the n-dimensional nor- 
mal distribution A/'„(zi*^* :2;^^t,t ,a:'j ^j^j^^t.t ,x ._ ^^i;t,t '^'j^^^^^ g^jid (j*'* .^: = 
{a ■" ' )l<ij<n- 
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For practical applications it is useful to discretize estimator (|23p to 



M 






</.(t,^-l).,i+,5(i,x-l).,i+,„0)-^^ ,*+..(*,.-...*+, ™«^y'^ i 

where [)i :— h{dii, . . . , (5j;„) ((5ij being the Kronecker symbol), for small enough 
h > 0. Without further details we note that according to Milstein and Tretyakov 
(2004) in a related context, it is efficient to take h « x/vM. 

As an alternative, it is also possible to expand the derivatives in (j23p , which 
leads to a SDE system of first order variation as in J28[ and [15| . In the differen- 
tiation of (|23p with respect to x for a fixed trajectory, r* '^(i^) can be considered 
to be independent of x. This can be seen as follows: if r* '^(w) = p, the random 
variable X* '^, which is assumed to have a density in M", lays almost surely in 
the interior of the exercise region. Due to the fact that (almost surely) the map 
X -^ Xp '^ is smooth (e.g. see Protter (1990)), X* '^ lays in the exercise region 
for y in an open disc around x. As a consequence, for any y in this disc we have 
T* '^{lo) — T, '^{lo) — p. Thus, by differentiating (|23p path-wise we obtain 



A 






?)■ 



1 Y^ P(^, X, ^+, g(^, X, ^+,„i 0) df t+.g(t,x,t+. „£;m 



dg{t,X,t+, rnCl 

T, ---- • — ' dx'- 



d p{t,x,t+,y) ag{t,x,t+, „;) 0/ 



where ^ ^(t'3:'t+ ) ' ' dx^' *" ' al '-^^^ ^^ principle be expressed analytically, 

and the vector process dyX^ '^(•) := — ^ — (•), s > t+, can in principle be 
simulated via a variational system of SDEs (e.g. see Protter [33|, Milstein & 
Schoenmakers [23], Giles & Glasserman 'iSl). 

In this paper we will prefer the discretized version (|26p of (|23p for our appli- 
cations. The algorithm is as follows. We first choose an /i > 0, and sample mS. 
for in = 1, . . . , M from the reference (usually normal) density. Next we simulate 
for each m a pair of trajectories mX^, which start in ^5^ := g{t,x±h, t~^, mO 

at i+, and end at the optimal stopping times m'''* '■= t* '™^ • Of course the 
optimal exercise dates ^t^ are generally unknown in practice, but we assume 
that we have good approximations mT at hand, which are constructed via 
some well known procedure. For example, in a pre-computation we may con- 
struct an exercise boundary via a regression method (e.g. LongstafF & Schwartz 
[27|), or as an alternative, we may use the policy iteration method of Kolodko 
& Schoenmakers [2J], see also Bender & Schoenmakers [5[- As discussed above, 
for a particular w we have mT^ — mT^ provided that h is small enough. For this 
reason we take in our simulations simply rnT~ =m t'^ , where r is some approx- 
imation of the optimal exercise policy. This pragmatic assumption is justified if 
the probability of the event mT~ ^^ r+ is small enough, i.e. h is small enough. 
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For each m we compute also the values mP^ '■= p{t, a; ± f), t'^, m5^) and m<j)^ 
:= (p{t, a; ± f), t~^ , my), and finally compute the estimate 



Remark 9. In the previous sections vector and matrix components are denoted 
by superscripts, so that time parameters of processes can be denoted by sub- 
scripts. In the next sections we depart from this convention and use subscripts 
for vector and matrix components. 

5 WKB approximations for transition densities 

5.1 Recap of WKB theory 

We summarize some results concerning WKB-expansions of parabolic equations 



(cf. Kampen [2l| for details). Let us consider the parabolic diffusion operator 

t + i«^t + 5E,,«..afe7 + S.&ii7- (28) 

For simplicity of notation and without loss of generality it is assumed that the 
diffusion coefficients a.y and the first order coefficients bi in ([28|l depend on the 
spatial variable x only. In the following let St := T — t, and let the functions 

{x, y) -^ d{x, y) > 0, [x, y) —> Ck{x, y), k>Q, 

be defined on K" x M", with (P and Ck, k > 0, being smooth. Then a set of 
(simplified) conditions sufficient for pointwise valid WKB-representations of the 
form 

p{t,x,T,y)^-^cxp(-^^^ + f2ck{x,y)6A, (29) 

V2^t \ 2St ^^ J 

for the solution (t, x) — > p{t, x, T, y) of the final value problem 

dp 

— — \- Lp — 0, with final value (30) 

piT,x,T,y) = <5(x-y), y G R" fixed, 
is given by 

(A) The operator L is uniformly elliptic in R", i.e. as in ^ the matrix norm 
of {aij{x)) is bounded below and above byO<A<A<oo uniformly in x, 

(B) the smooth functions x — > aij (x) and x —f bi (x) and all their derivatives 
are bounded. 



For more subtle (and partially weaker conditions) we refer to Kampen 2lj. If 
we add the uniform boundedness condition 

(C) there exists a constant c such that for each multiindex a and for all 1 < 
i,j,k < n, 



d°'a 



jk , 



dx° 



d"b. 



dx° 



(^)' 7^(^) <cexp(c|xp), (31) 
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then the Taylor expansions of the functions d and Ck around y G M" are equal 
to d and Ck,k > globally. I.e. we have the power series representations 

d^{x,y) = ^d„(2/)fe« (32) 

Ck{x,y) = ^Ck,a{y)Sx°', k>0, (33) 

a 

where 6x := x — y. Note that (C) is implied by the stronger condition that all 
derivatives in (j3ip have a uniform bound. Summing up we have the following 
theorem: 

Theorem 10. // the hypotheses (A), (B) are satisfied, then the fundamental 
solution p has the representation 

p{St,x,y) = -^exp l~^^^ + Y,Ck{x,y)StA , (34) 

where d and Ck are smooth functions, which are unique global solutions of the 
first order differential equations p5p .()36 p . and psp below. Especially, 

i5t,x,y) ^ Stlnp{St,x,y) = ~^Stliii2TrSt) - ^ +J2ck{x,y)St''+^ 

fc>0 

.2 

is a smooth function which converges to —-^ as St ^ 0, where d is the Rieman- 
nian distance induced by the line element ds^ = ^^^ a~j dxidxj, where with a 

slight abuse of notation (a~- ) denotes the matrix inverse of{aij). If the hypothe- 
ses (A),(B) and (C) are satisfied, then in addition the functions d, Cfc,fc > 
equal their Taylor expansion around y globally, i.e. we have i32\) - i3S\) . 

The recursion formulas for d and Ck, k > are obtained by plugging the 
ansatz (|29p into the parabolic equation ((30|) . and ordering terms with respect 
to the monoms St'' = (T — <)* for i > —2. By collecting terms of order St^^ we 
obtain 



"' 4 



J^'^l^'^ijdl^, (35) 



where d^ denotes the derivative of the function d^ with respect to the variable 
Xk, with the boundary condition d[x, y) = for x = y. Collecting terms of order 
Sf^^ yields 

- f + ^Ld' + ^E IE K(^) +«^'(^)) %) |^(^'2^) = 0' (36) 

where the boundary condition 

co(y,y) = -2ln\/det(ay(y)) (37) 

determines cq uniquely for each y e M". Finally, for fc + 1 > 1 we obtain 

{k + l)ck+^{x,y) + i^^^.a,,(x)(^%±i + %%f ) 

(38) 
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with boundary conditions 

Ck+i{x,y) = Rk{y,y) a x = y, (39) 

Rk being the right side of p8|) . For some classical models in finance a global 
transformation of the diffusion operator to the Laplace operator is possible (at 
the price of more complicated first order terms however). We observe this in 
the case of a Libor market model (Section E]). The requirement that a transfor- 
mation y{x) of an operator with second order coefficients aij{x) = {aa^)ij{x) 
leads to a Laplacian with respect to second order terms in y is equivalent to 

^(aa Ui,) — ^^^S,k (40) 

ml 

where Sjk denotes the Kronecker delta. If cr is invertible it follows directly that 
the transformation y{x) satisfies the first order matrix equation 

(^) = K7(-))- (41) 

The latter equation determines the transformation (up to constants, of course) 
but cannot be integrated in general, and if it can not explicitly in general. 
However, a necessary and sufficient condition for integrability of ()4ip in terms 
of a can be given, where we restrict ourselves to the case of invertible cr. 

Proposition 11. There is a global coordinate transformation for the opera- 
tor I128\} such that the second order part of the transformed operator equals the 
Laplacian, iff Oij = {aa^)ij for a (square) matrix function a which satisfies 

-—^ai,ix)^}^-^aikix), xeR\ (42) 

1=1 ' 1=1 ' 

The latter fact is also observed and proved in Ait-Sahalia [1|. If the condition 
of Proposition [TT] is satisfied, then coordinate transformation leads to second 
order coefficients of the form a^ = Sij , so that the solution of ([55]) becomes 

d'^{x,y) = ^{xi -yif- 

i 

If conditions (A), (B), (C), and (|^^ hold, then in the transformed coordinates, 
explicit formulas for the coefficient functions Ck,k > can be computed via the 
formulas 

Co{x,y) = ^{yi-Xi) I bi{y + s{x~y))ds, 

Jo 

Ck+i{x,y) = / Rk{y + s{x-y),y)s''ds, (43) 

Jo 

with Rk being the right-hand-side of (J38D where a.y = Sij . Similar formulas are 
obtained in Ait-Sahalia Jj. In Kampen [21i] it is shown in addition how the co- 
efficients Ck can be computed explicitly in terms of power series approximations 
of the diffusion coefficients Oij and 5,. However, in high dimensional models 
such as the Libor market model direct computation of the coefficients Ck seems 
more feasible as it turns out that the computation up to the coefficient ci is 
sufficient for our purposes. 
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5.2 Error estimates 

We now study the approximation error of a truncated WKB expansion (and 
its derivatives), which is essential for convergence of the Monte Carlo schemes. 
In this respect we will show how the derivatives (up to second order) of the 
product value function with respect to the underlyings computed by means of a 
truncated WKB-expansion converge in supremum norm and Holder norms. Let 
us consider a WKB-approximation of the fundamental solution p of the form 

pi{t, X, T, y) = — i=^ exp I -^-^ + E ^^i^' y)^*'] ' (44) 

i.e. we assume that the coefficients (P and Ck, < k < I have been computed 
up to order I (recall that St — T — t ioi the sake of brevity). Let us denote the 
domain of the Cauchy problem by I? = (0, T) x M". For integers n > and real 
numbers S G (0, 1) let (7™+'5/2,n+<5(|£)^ \)q ^^q space of m (n) times differentiable 
functions such that the mth (nth) derivative with respect to time (space) is 
Holder continuous with exponent ^ (S). Furthermore, |.|m+i5/2.n+(5 denote the 
natural norms associated with these function spaces. It is well-known that in 
case of our assumptions (A) and (B) the fundamental solution p satisfies the a 
priori estimate 

\p{t,x,T,y)\ < C{T-tr^^/'e^p (- ^2(r-tf ) ' ^^^^ 

for some generic constant C and some Aq which is less or equal than the lower 
bound A in assumption (A) above. We call a WKB-approximation pi of the 
fundamental solution p admissible, if it satisfies the a priori estimate (j45p . The 
WKB-approximation po is always admissible while the proof of theorem 9 (cf. 
[2l|) that pi is admissible if / > ^o where Iq is some natural number depending 
on the coefficient functions and can be computed by comparison of the WKB- 
expansion and the Levy-expansion. For lower I admissibility has to be ensured 
for each model. In the Libor market model admissibility for / = 1 is ensured. 
As a consequence of Safanov's theorem (cf. Krylov [25|) we have 

Theorem 12. Assume that (A), (B), and (C) are satisfied and let h G C^+'^ (R") 
and f e C^/^^\D). If 

c < —A for some A > 0, (46) 

then the Cauchy problem 



W + hJ:^,avi^)ol7Bt+J:^b^ix)§^^+cix)w = f{t,x) zn D 



(47) 



w{T, x) = h{x) for X G 



has a unique solution w, and there exists a constant c depending only on S, n 
A, A anrf /iT = max{|a|5, |6|a-, jcj^} such that 

\w\l+S/2,2+5 < c [\f\s/2.s + \h\2+s] ■ (48) 

In order to analyze the truncation error of the Cauchy problem with data h 
we consider the function 

u {t,x) — u{t,x) ~ ui{t,x), where 
15 



u{t,x)= h{y)p{t,x,T,y) and ui{t,x)^ h{y)pi{t,x,T,y)dy. (49) 

We say that ui in ([1^ is admissible if p; is admissible. It is now possible to 
derive different error estimates in strong norms depending on which Greeks we 
want to control on which level of regularity. 

Theorem 13. Assume that conditions (A), (B), and (C) hold and that h G 
C'^^ (M") and assume that ui is admissible. Then 

\u{t,x) - ui{t,x)\i+s/2,2+s e 0(t'"^). 

Proof. Let w{t,x) — e^^^u^{t,x) with r constant and wi{t,x) — e^^*ui{t,x). 
Since 

^+Lu^ = -^-Lui =: f^, (t, x), we have 

^+Lw + rw = ^ + Lwi+rwi^ e^* ("ff " ^ui) =: e'-*/„, (i, x). 

Admissibility of ui and an argument similar to that of Krylov ensures that the 
right side of the latter equation can be measured in the norm |.|5/2,5- Hence we 
can apply the estimate (|48p to the function w{t, x) = e~^*u^{t, x) for a constant 
r > and we get (after dividing by e''*) 

l"'^ 1 1+5/2,2+5 < C\fui\s/2,S- (50) 

In order to compute the term on the right side of ([50]) we can plug (|ii| into the 
left-hand side of PUJ) the parabolic equation satisfied by the exact fundamental 
solution p. However in order to see how the higher order terms behave exactly 
we plug in 

1 / d^(x v) ' \ 

p{t, x, y) = — =^ exp yi^ + J2 Cfc(a^, y)t'" + Ri+i{t. x,y)\, (51) 



where 



We get 



Ri+i{t,x,y) = Y. ^k{x,y)t'' = 0(i'+i). (52) 

k=l+l 



at ^ 2 Z^y "-V aa;i9xj ^ l^i "* dxi ~ 

+(q,x, + i?i+i,^J(Q,^^. + i?/+i,x,)) + i(Q + Ri)\p = 0{t')p, t i 0. 

Applying a priori estimates for p we get the result. ■ 

Remark 14. A more intricate analysis shows that for practical purposes it 
is possible to remove the admissibility condition above, if we approximate the 
Cauchy problem by a Dirichlet problem with a large but spatially bounded 
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domain (a natural step from the numerical point of view). However, since 
this involves an additional analysis of an integral equation corresponding to the 
boundary condition we go not into the details here. Generalizations to estimates 
which include Taylor-expansions of the WKB-coefRcients will also be considered 
elsewhere. The assumption h S C^+''(R") can be weakened to Holder continuous 
pay-offs if we abstain from controlling the Q Greek (sensitivity with respect to 
time) up to maturity. The case where Ck are computed up to fc = 1 is the first 
case where the truncation error for first and second derivatives converges to zero 
(in supremum norm with order 0{St) and in Holder- extension of supremum 
norm with order 0{Sty^^). This implies that our Monte Carlo computation 
scheme for the Greeks converges. 

Remark 15. We can easily see how the boundedness of the constant M5 in 
Theorem|4]is controled in situations where p is WKB-approximated by p', ^ > 0, 
and where the prior is chosen according to (|24p in Section 4. For simplicity 
we here assume that the problem is reduced to the form a^j = 6ij. Then the 
logarithmic derivative of the WKB-expansion (cf. (35)) of the density p takes 
the form 

The logarithmic derivative of the lognormal prior (|24p takes the form 

i.e., in the difference the first terms cancel out. So for small time t the main 
contibution to the constant M5 is the difference of ^(cg — %) which does not 
depend on t. Note that if we freeze (|53p at xq say, which leads to the naive 
estimator ([TTj) . the difference contains a term of order of 0{t~^)\ 

6 Applications to the Libor market model 

We consider a Libor market model with respect to a tenor structure < Ti . . . < 
Tn+i in the terminal measure Pn+i (induced by the terminal zero coupon bond 
Bn+i{t)). The dynamics of the forward Libors Li{t), defined in the interval 
[0,Ti] for 1 < i < n, are governed by the following system of SDE's (e.g., see 
Jamshidian [20|), 






j=i+l ■> ■' 

(54) 
where 6i — Ti+i—Ti are day count fractions and t -^ Ji{t) = (ji.iit), ■ • ■ , 7i,d(0); 
<t < Ti, are bounded and smooth enough deterministic volatility vector func- 
tions. We denote the matrix with rows 'yj by F and assume that F is invertible. 
In what follows we assume that F(t) = F does not depend on t. The case of 
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time-dependent volatility is discussed in Kampen, Kolodko, Schoenmakers |22| . 
In ((5l)) . (iy("+^^(t) I < i < r„) is a standard d-dimensional Wiener process 
under the measure P„+i with d, 1 < d < n, being the number of driving factors. 
In what follows we consider the full-factor Libor model with rf = n in the time 
interval [0, Ti). 

6.1 WKB approximations for the Libor kernel 



Let us transform the dynamics of (|54p to Ki :— In Li, 1 < i <n, 

dK, = ^dU - l^d{U) = (-^ + A^»(^' e^S ■ • • , e^")) dt + ij dW^^+^^ , 

(55) 
where fi is given in ([5l|) . Note that the coefficients of the generator correspond- 
ing to process K are bounded and satisfy the WKB assumptions (A) (B) in 
Section 5.1. Hence, we may apply a WKB approximation to the transition 
density of the process ([55|l . 

By the transformation Y := T^^K we obtain the process 

dYi = ^f {t, Y)dt + dVK/"+^^ , 1 < i < n, where (56) 

n " I 12 



i=i i=i 

for which the generator has a Laplacian diffusion term, which leads to technically 
more convenient expressions in the respective WKB expansion. 

The situation of time independent 7 (hence bounded ^) is exemplified 
in case study Section 16. 2[ where the transition density p^ is approximated and 
subsequently transformed to an approximated transition density p^ of the Libor 
process. Below we spell out the ingredients for computing the corresponding 
WKB coefficients cq and ci according to (|15)) to be exploited in Section 16.21 
Using the notations 

^'^^'"'^)^= (r(.-,)), ^" l + ^.e(r.V ^^^^"' 
and a :— {'yj'yj)2j^ij we may write, 

n n 71 n 

co{s,x,y) =^Vi{y,-Xi) + ^^Tr.^{y,-Xi) ^ ajiFi{s,x,y), (57) 

i=l i=l j=l l=j+l 

o n n 



dxp 



j=i i=j+i 

dFi{s,x,y) 



d^co, . oV^T^-i V^ dFi{s,x,y) 



^ j^l /^j + 1 ^ 

n n n 



Y^Y^T.-!. N Y^ d'^Fi{s,x,y) 

i=i j=i i=j+i P 
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where 



dFi{s,x,y) 
dxr, 



r;p(s) 



{T{x-y))i \1 + Sie(^^)' 



<5,e(r-) 



Fiis,x,y), 



and 



d^Fi{s,x,y) 
dx'i 



FAsxy)^^!^:^ 



2r?p 






S,e(^-h 



{T{x-y))i{l + Sie(^-)'r- 



We finally obtain p^{s, u, t, v) by density transformation formula, 



p''(s,u,t,v)^p'^{s,Sz\u),t,Sr\v)) 



dSi\v) 



dv 



with 



Si\v):=T-Ht){lnvi,...Mvny 



For simplicity in the case study below we assume that the matrix F is upper 
triangular and does not depend on t. We then have, 



p^{s,u,t,v) 



VMt - s) 



" P-I / F-i(ln^,...,ln^)^ 



2(t-s) 



-Y,Ckis,S-'{u),S-'{v))it-s)'^ 



k=0 



6.2 Case study 

We now illustrate the estimators ([9]) and ([T5|) in Section [3] and the estimators 
(PIj) and (P5|) in Section [J] by computing European and Bermudan swaptions 
and Deltas in a Libor market model. A (payer) swaption contract with maturity 
Ti and strike with principal $1 gives the right to contract at Tj for paying a 
fixed coupon 9 and receiving floating Libor at the settlement dates T^+i,. . . ,r„. 
The discounted payofi of the contract is thus given by 



MLiT,)) 



Bn+l{Ti) ^ 



Y,Bj+i{T,){S,L,{T.,) 



(58) 



j=i 



For our experiments we take in ([5l)) . Si = 0.5, L{0) — 3.5% flat, and constant 
volatility loadings, 7i(t) = 0.2ei, where e^ are n-dimensional unit vectors de- 
composing an input correlation matrix p. 



Pij = exp 



n- 1 



Inpo 



1 < j,j < n. 



(59) 



with Poo = 0.3 (for more general correlation structures we refer to Schoenmakers 
[35J). We consider at-the-mondey European swaptions with maturity Ti and at- 
the-money Bermudan swaptions with 10 annual exercise possibilities, starting 
from Ti, hence 9 — 3.5% in 
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For the Bermudan swaptions a good approximation of the optimal stopping 
pohcy is constructed by Andersen's method (strategy II, see Andersen ^), 



Ti,L(0) 



inf {i : T,>Ti,T,e T, B„+i(T,)/^(4j''^'"') > H,+ 



B„+i(r,)/,(4;'^(°)) \, (60) 



max E^ 

j>i,T,e T 

where L°j.^^°^'^ := Lj{T^) in hnc with Sections [D and H and T := {Ti, T3, T5, 
Tig} is the set of possible exercise dates. The conditional expectations in 



((50)) can be computed accurately in closed- form (see, e.g., Schoenmakers j35l|). 
Further in (|60p. H is a constant vector computed by backward optimization over 
a set of pre-simulated trajectories, as proposed by Andersen 2]. In Table 3, 
column 2, we display the Bermudan prices u|,°" due to stopping strategy ta- 
Upper estimations u"^ are constructed from v}°^ by the dual approach, developed 
in Rogers 34| and Haugh and Kogan [1»|, see Table 3, column 1. As we see, the 



distance between lower and upper Bermudan estimates does not exceed 0.5% 
(relative to the values). 

The Libor transition kernel p^(s, x, t, y) shows to have a pronounced "delta- 
shaped" form. Because of this, it is very important for efficiency of the esti- 
mators in Sections [3]|1] to find a suitable proxy density (j). We take for <j) the 
transition kernel of a lognormal approximation L''^", obtained from the Libor 
process ()54p by freezing the coefficients at the initial time s. 



Lt "'""''(O =2:iexp(^,), (61) 

where ^ is a n-dimensional Gaussian vector with 



Co«(C.,0) = r,,, l<i,3<n. (62) 



The transition density of L'"" is then given by 



1 -rV rr^ 






exp 



r-i((lnf^...ln^)-Ai'-(s,i,x))^ 



2{t~ s) 



with I • I denoting the Euclidean norm. So, in order to sample from density 0, 
we simulate via (pT |) - (|5^ the lognormal samples 

,„C = 4"°''^'"\mO=:5(0,i(0),ri,„O, m = l,...,M. 

As a (more accurately approximated) Libor transition kernel, we use WKB 
approximation pg and pf . We endow the corresponding estimators with super- 
scripts and 1 respectively. 

European and Bermudan prices and Deltas via the estimators © , ([TS]) , (PTjl , 
(|26p are given in Tables 1-4. These results are compared with corresponding 
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estimates due to "exact" Libor trajectories, simulated by a log-Euler scheme 
with smah time step At (we take At = 5i/b for Europeans, At = 5j/10 for 
Bermudans). The "exact" estimates are endowed with the superscript ox. For 
comparison, the corresponding estimates due to the standard lognormal Libor 
approximation L'"*" are computed as welL In order to keep standard deviations 
within 0.5% relative (to the values) we take /i = 3.5 x 10~^, M = 5 x 10^. As we 
see, the WKB approximation with only two coefficients, cg and ci, provides a 
very close estimate of the European swaptions and Deltas, also for large matu- 
rities. The distance between the values simulated via "exact" Libor trajectories 
and the corresponding values due to the WKB approximation is smaller than 

0.5% relative to the value. In contrast, the lognormal estimators /ig„, g''^" , Uig„ 



and 



dxi 



give an acceptable approximation only for Ti < 2. 



Table 1. European 


swaptions (values in bas 


s points) 


Ti 


/e. (SD) 


/,,„ (SD) 


h (SD) 


h (SD) 


1.0 
2.0 
5.0 
10.0 


178.9(0.4) 
245.3(0.6) 
351.3(1.0) 
429.6(1.5) 


179.0(0.4) 
246.5(0.6) 
359.8(1.0) 
451.4(1.6) 


181.6(0.4) 
251.4(0.6) 
376.4(1.1) 
495.6(1.7) 


178.9(0.4) 
244.3(0.6) 
352.7(1.0) 
431.8(1.4) 



Table 2. European Deltas (values in 


basis points) 




Ti 


fr (SD) 




- — -Ci) 

it (SD) 


§Z (SD) 


1.0 
2.0 
5.0 
10.0 


1768.3(2.8) 
1726.4(2.9) 
1599.6(3.2) 
1417.1(3.8) 


1774.2(2.8) 
1732.1(2.9) 
1615.9(3.3) 
1474.0(4.2) 


1794.9(2.9) 
1729.0(2.9) 
1722.5(3.5) 
1668.0(4.7) 


1770.7(2.8) 
1729.0(2.9) 
1597.0(3.2) 
1422.7(3.9) 



Table 3. Bcrmudan swaptions (values in basis points) 



Ti 


w'7 (SD) 


ui; (SD) 


K. (SD) 


uo (SD) 


ui (SD) 


1.0 


351.2(0.7) 


352.5(1.0) 


350.9(0.7) 


354.7(0.7) 


351.2(0.7) 


2.0 


388.4(0.8) 


389.8(1.0) 


388.2(0.8) 


396.6(0.8) 


387.3(0.8) 


5.0 


461.5(1.1) 


463.4(1.3) 


466.3(1.1) 


492.9(1.1) 


460.8(1.1) 


10.0 


523.7(1.6) 


524.8(1.7) 


543.6(1.7) 


601.2(1.7) 


523.6(1.5) 



Table 4. Bcrmudan Deltas (values in 


basis points) 




Ti 


tf (SD) 


^ (SD) 


tr (SD) 


ir (SD) 


1.0 
2.0 
5.0 
10.0 


2709.2(3.5) 
2631.1(3.5) 
2392.9(3.7) 
2101.5(4.4) 


2720.9(3.5) 
2630.5(3.5) 

2407.7(3.8) 
2152.5(4.7) 


2747.2(3.5) 
2700.7(3.6) 
2561.9(4.0) 
2443.4(5.3) 


2709.2(3.5) 
2628.6(3.5) 
2398.0(3.8) 
2111.5(4.4) 



Remark 16. The values in Tables 1-4 are computed using a second order Taylor 
approximation of ci (x, y) around x, where ci (x, x) , the derivatives ^^ (x, x) and 

d 'd (^;^)i ^^6 computed (using finite differences) prior to the Monte Carlo 
simulation. 
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Figure 1: CPU time (seconds) for simulating European (left) and Bermudan 
(right) Deltas for different Ti by log-Euler Libor simulation (solid line) and by 
WKB density approximation (with cq and ci) (dash line). 

Computational time 

By using the new estimators we avoid step-by-step Euler simulation of the Libor 
process in the time interval [0, Ti]. Generally, the cost of Euler stepping up to Ti 
is proportional to Ti/At, whereas the cost of the "direct estimators" P^ and 
(j26p is independent of Ti. In particular, in the present Libor case, Euler stepping 
up to Ti requires a cost proportional to n^ -^ times the cost of computing the 
(possibly virtual) pay-off at Ti. In comparison, the cost of simulating estimators 
(fT5|l and ((26|) is proportional to ri^ times the cost of the pay-off at Ti . 

In Figure 2 we compare for different Ti the CPU time (per sample) needed 
for computing the values in Tables 2,4 using WKB based estimators (|15p and 
([2^ with the CPU time required for computing the estimates via straightfor- 
ward Euler stepping of "exact" Libor trajectories up to Ti. We conclude that, 
particularly for larger Ti, the efficiency gain is quite high in the European case, 
and still considerable in the Bermudan case. 

Remark 17. Fries & Kampen 14 1 and Fries & Joshi [13'| propose simulation 
schemes which improve upon Euler SDE simulation and allow for taking larger 
time steps for obtaining the same accuracy. Assuming that such a scheme re- 
quires a time step of order, say 0(-\/At), instead of 0{At) for the same accuracy, 
it is clear that, for example in the European case, the gain of our method with 
respect to this one is still order of 0(T/^/At). 
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